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Abstract 

A  new  finite  element  method  is  presented  that  features  the  ability  to  include 
in  the  finite  element  space  knowledge  about  the  partial  differential  equation 
being  solved.  This  new  method  can  therefore  be  more  efficient  than  the  usual 
finite  element  methods.  An  additional  feature  of  the  partition-of-unity  finite 
element  method  is  that  finite  element  spaces  of  any  desired  regularity  can  be 
constructed  very  easily.  Moreover,  the  method  is  of  “meshless”  type.  This 
paper  includes  a  convergence  proof  of  this  method  and  illustrates  its  efficiency 
by  an  application  to  the  Helmholtz  equation  for  high  wave  numbers.  The  basic 
estimates  for  a-posteriori  error  estimation  for  this  new  method  are  also  proved. 
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1  Introduction 


In  this  paper,  we  present  a  new  method,  the  “partition  of  unity  finite  element  method” 
(PUFEM)  to  construct  conforming  finite  element  spaces  (FE  spaces)  with  local  prop¬ 
erties  determined  by  the  user.  The  features  of  the  PUFEM  are: 

1.  The  ability  to  include  a  -priori  knowledge  about  the  partial  differential  equation 
in  the  FE  space, 

2.  the  construction  of  FE  spaces  of  any  given  regularity  (which  is  a  desirable  feature 
for  the  approximation  of  higher  order  equations), 

3.  the  method  is  “mesh-free”, 

4.  the  method  can  be  understood  as  a  generalization  of  the  classical  finite  element 
methods;  in  particular  the  h,  p,  and  hp  versions  of  the  finite  element  method 
can  be  understood  as  special  cases  of  the  PUFEM, 

5.  the  PUFEM  permits  a  -posteriori  error  estimation  and  adaptive  approaches. 

Let  us  elaborate  these  five  features  in  more  detail.  The  first  feature  touches  the  ques¬ 
tion  of  approximation  properties  of  the  FE  spaces.  The  classical  FE  spaces,  known 
as  h  and  p  version  of  the  finite  element  method,  are  spaces  which  have  good  lo¬ 
cal  approximation  properties  and  are  conforming;  typically  they  consist  of  piecewise 
polynomials  (or  mapped  polynomials)  and  satisfy  some  continuity  requirement  across 
inter-element  boundaries.  In  the  h  version,  the  polynomial  degree  is  fixed  (typically 
p  <  2)  and  approximation  is  achieved  by  decreasing  the  meshsize  h.  An  appropriate 
interpolant  (e.g.,  for  p  =  1  on  triangles,  nodal  interpolation  can  be  taken)  produces 
a  good  approximation  which  satisfies  the  necessary  continuity  conditions.  In  the  p 
version,  local  approximation  is  realized  by  polynomials  of  increasingly  higher  degree. 
The  approximation  properties  of  conforming  p  extensions  are  due  to  two  facts.  Un¬ 
constrained,  i.e.,  without  any  inter-element  continuity  constraints,  polynomials  have 
good  approximation  properties  on  each  patch.  The  resulting  jumps  across  inter¬ 
element  boundaries  can  be  resolved  by  polynomial  corrections  because  polynomial 
spaces  are  large  enough  to  permit  continuous  extensions  from  the  element  boundaries 
into  the  elements  (see  [18], [20]). 

For  many  problems,  the  particular  structure  of  the  equation  can  be  exploited  to 
construct  local  spaces  with  better  approximation  properties  than  the  usual  h  or  p 
version  based  spaces.  For  example,  solutions  to  Laplace’s  equation  (A u  =  0)  in  the 
plane  can  be  approximated  by  harmonic  polynomials  only;  it  is  not  necessary  to  use 
all  polynomials  for  an  approximation  in  a  p  version  fashion.  However,  as  opposed 
to  full  polynomial  spaces,  there  are  not  enough  harmonic  polynomials  to  construct 
conforming  spaces  which  consist  of  piecewise  harmonic  polynomials.  The  PUFEM 
however,  allows  us  to  construct  conforming  spaces  from  harmonic  polynomials  and 
thus  exploit  their  good  approximation  properties. 


Let  us  mention  a  few  examples,  for  which  the  special  structure  of  the  underlying 
partial  differential  equation  can  be  exploited.  “Generalized  harmonic  polynomials” 
can  be  constructed  for  elliptic  equation  with  analytic  coefficients.  These  “general¬ 
ized  harmonic  polynomials”  have  approximation  properties  very  similar  to  those  of 
harmonic  polynomials  (see  [23,  3,  7,  8,  11]).  For  Laplace’s  equation  or  the  elasticity 
equations,  corners  or  sudden  changes  of  boundary  conditions  produce  certain  types  of 
singularities  that  are  resolved  very  poorly  by  the  classical  methods  unless  a  properly 
refined  mesh  is  used.  [14,  15]  show  that  the  use  of  special  shape  functions  for  these 
problems  can  be  very  efficient. 

For  problems  of  3-d  elasticity  on  polyhedral  domains,  the  creation  of  a  properly 
refined  mesh  can  be  complicated  in  the  vertices,  where  vertex  and  several  edge  singu¬ 
larities  interact.  The  PUFEM  circumvents  this  difficulty  because  it  enables  the  user 
to  incorporate  the  singularities  directly  in  the  FE  space. 

The  classical  methods  work  well  if  the  solution  to  be  approximated  is  smooth.  How¬ 
ever,  for  many  problems  of  practical  interest,  the  solution  is  not  smooth.  It  can  be 
highly  oscillatory  as  in  the  case  of  Helmholtz’s  equation  with  high  wave  numbers, 
or  it  can  be  rough,  as  in  the  case  of  the  analysis  of  composites,  laminated  materials 
or  stiffeners.  The  usual  finite  element  methods  can  be  prohibitively  expense  for  this 
kind  of  problems.  However,  as  shown  in  [17],  the  use  of  shape  functions  reflecting 
this  rough  behavior  can  lead  to  optimal  convergence  rates.  A  similar  observation 
can  be  made  for  highly  oscillatory  problems  such  as  Helmholtz’s  equation.  It  was 
demonstrated  in  [11]  that  the  approximation  with  plane  waves  displaying  the  same 
oscillatory  behavior  as  the  solution  is  very  efficient. 

Another  example,  where  non-polynomial  approximation  spaces  are  of  physical  inter¬ 
est,  are  problems  on  unbounded  domains.  For  problems  such  as  Laplace’s  equation 
or  Helmholtz’s  equation,  expansions  of  the  solution  around  the  point  at  infinity  are 
known  and  can  be  used. 

Let  us  comment  on  the  second  feature  of  the  method,  the  ability  to  construct  FE 
spaces  of  any  given  regularity.  The  PUFEM  creates  FE  spaces  as  follows  (a  more 
detailed  description  follows  below).  Let  patches  {fl,}  comprise  a  covering  of  the 
domain  f2,  and  let  {95,}  be  a  partition  of  unity  subordinate  to  the  covering.  On  each 
patch,  let  function  spaces  Vi  reflect  the  local  approximability.  Then  the  global  finite 
element  space  V  is  given  by  V  =  YT  (piVi.  Local  approximation  in  the  spaces  V.  can 
either  be  achieved  by  the  smallness  of  the  patch  (an  h  version)  or  by  good  properties  of 
Vi  (a  p  version).  Theorem  1  states  that  the  global  space  V  inherits  the  approximation 
properties  of  the  local  spaces  Vi.  Addionally,  it  inherits  the  smoothness  of  the  partition 
of  unity  (and  the  spaces  V^).  Therefore,  the  construction  of  smoother  FE  spaces  for  the 
approximation  of  higher  order  equations  as  they  appear,  for  example,  in  various  plate 
and  shell  models  is  easily  possible  by  using  a  partition  of  unity  which  is  sufficiently 
smooth. 

Let  us  turn  to  the  third  feature.  The  method  is  mesh-free  in  the  sense  that  no  mesh 
has  to  be  generated  explicitly;  rather,  a  “mesh”  is  determined  implicitly  through  the 
overlapping  patches  of  the  covering.  Let  us  observe  that  the  creation  of  a  partition 
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of  unity  can  be  easy  (if  W{  are  functions  living  on  the  patches  fi,,  the  functions 
<Pi  =  Wil  Yhjwi  f°rm  a  partition  of  unity)  and  fully  automated.  Hence,  changing  the 
FE  space  adaptively  is  easy.  Moreover,  changes  as  they  are  very  common  for  problems 
like  crack  propagation  or  the  optimal  placement  of  a  fastener  can  be  done  easily.  In 
the  latter  example,  for  instance,  the  user  would  like  to  calculate  several  possible 
locations  of  the  fastener.  In  the  usual  finite  element  method,  he  has  to  remesh  (at 
least  locally)  for  each  case,  which  could  be  costly.  In  the  PUFEM,  the  effect  of  the 
fastener  could  be  modeled  by  a  special  function  and  hence,  for  each  run,  only  a  few 
spaces  Vi  (in  which  these  special  function  are  contained)  have  to  be  changed.  “Mesh- 
free”  methods  have  been  proposed  recently  in  [9,  13].  For  an  appropriate  choice  of 
parameter  values,  the  method  of  [9]  reduces  to  a  particular  type  of  PUFEM,  and  the 
computational  analysis  of  [13]  shows  that  the  method  of  [9]  is  most  efficient  in  that 
case. 

Commenting  on  the  fourth  point  of  the  above  list,  we  observe  that  the  if  we  ap¬ 
proximate  locally  on  patches  by  polynomials  of  degree  p,  we  get  a  FE  space  with 
approximation  properties  similar  to  the  usual  finite  element  methods.  If  the  approx¬ 
imation  in  the  spaces  U  is  achieved  through  the  smallness  of  the  patches,  we  get  a 
method  very  similar  to  the  h  version;  if  approximability  is  realized  by  an  increase  of 
the  polynomial  degree,  the  method  behaves  like  the  p  version.  If  both  are  varied,  we 
get  an  hp  version.  In  this  sense,  the  method  presented  here  is  a  generalization  of  the 
usual  finite  element  method.  It  has  the  feature  to  include  in  the  FE  space  knowledge 
about  the  structure  of  the  particular  problem  at  hand;  however,  one  can  achieve  local 
approximation  by  polynomials  and  then  the  method  produces  essentially  the  usual 
FE  spaces. 

Concerning  the  fifth  point  of  the  list,  we  mention  that  the  PUFEM  permits  a- 
posteriori  error  estimation.  The  basic  ideas  for  a  -posteriori  error  estimation  were 
developed  in  [1],  [2]  and  in  our  development  of  a-posteriori  error  estimation,  we  will 
follow  closely  [1]. 

For  a  successful  implementation  of  the  PUFEM,  three  issues  have  to  be  addressed: 

1.  The  integration  of  the  shape  functions  constructed  by  the  PUFEM. 

2.  Finding  a  basis  of  the  PUFEM  space  and  controlling  the  condition  number  of 
the  stiffness  matrix  created  by  the  PUFEM. 

3.  The  implementation  of  essential  boundary  conditions. 

Let  us  just  briefly  outline  why  these  three  issues  arise.  A  more  detailed  analysis 
of  these  issues  will  be  done  in  subsequent  work.  On  the  integration  issue,  we  ob¬ 
serve  that  typically  the  shape  functions  are  defined  on  the  patches.  Integrating  two 
shape  functions  against  each  other  requires  an  integration  over  the  intersection  of  two 
patches.  Hence  the  integrator  has  to  be  able  to  integrate  efficiently  over  intersections 
of  patches.  The  issue  of  finding  bases  of  the  PUFEM  spaces  and  the  problem  of  con¬ 
trolling  the  conditioning  numbers  of  the  stiffness  matrices  is  illustrated  in  section  3.3. 
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Finally,  let  us  state  that  essential  boundary  conditions  can  be  implemented  in  differ¬ 
ent  ways.  Either  the  local  approximation  spaces  on  patches  close  to  the  boundary 
are  chosen  large  enough  to  permit  an  extension  of  the  boundary  data  from  boundary 
into  the  domain,  or  Lagrange  multiplier  or  penalty  methods  are  used. 

The  paper  is  organized  as  follows.  In  section  2  we  develop  the  PUFEM  and  give  a 
proof  of  its  approximation  properties.  In  sections  3. 1-3.4,  we  illustrate  some  of  the 
features  of  the  PUFEM  in  a  one  dimensional  setting.  In  section  3.1,  we  demonstrate 
how  the  PUFEM  produces  robust  finite  element  spaces  for  a  problem  with  a  boundary 
layer.  The  performance  of  the  PUFEM  for  this  particular  problem  is  comparable 
to  the  usual  finite  element  methods  for  problems  with  smooth  solutions  because  the 
PUFEM  allows  us  to  create  finite  element  spaces  which  capture  precisely  the  behavior 
of  the  boundary  layer.  Section  3.2  proposes  several  types  of  partitions  of  unity  which 
satisfy  the  necessary  conditions  for  the  PUFEM  to  work.  Section  3.3  analyzes  in 
more  detail  the  case  of  polynomial  local  approximation  spaces.  In  particular,  the 
problem  of  potential  linear  dependencies  and  the  issue  of  the  condition  number  of 
the  stiffness  matrix  is  addressed.  In  section  3.4  finally,  a  PUFEM  is  exhibited,  in 
which  all  the  degrees  of  freedom  have  the  meaning  of  the  value  of  the  approximating 
function  in  appropriate  points.  Sections  4.1  and  4.2  discuss  briefly  methods  how  to 
choose  good  local  approximation  spaces  and  the  issue  of  the  optimality  of  local  spaces. 
Two  numerical  examples  are  presented.  In  section  6.1,  the  PUFEM  is  compared  with 
the  usual  p  versions  for  the  approximation  of  harmonic  functions.  In  section  6.2 
the  PUFEM  is  used  for  the  approximation  of  solutions  to  Helmholtz’s  equation  with 
large  wave  number.  The  PUFEM  is  shown  to  be  superior  (both  in  terms  of  error  per 
degree  of  freedom  and  error  per  floating  point  operation)  to  several  h  version  type 
finite  element  methods.  The  paper  concludes  with  a  proof  of  an  a  posteriori  error 
estimator  for  the  PUFEM,  which  is  based  on  exact  solutions  of  appropriate  local 
problems. 

2  The  Method 

In  this  section,  we  present  our  method  of  constructing  conforming  subspaces  of  If1  (SI). 
We  construct  FE  spaces  which  are  subspaces  of  ff1(S 2)  as  an  example  because  of  their 
importance  in  applications.  We  would  like  to  stress  again  that  the  method  leads  to 
the  construction  of  smoother  spaces  in  a  straight  forward  manner.  Crucial  to  the 
construction  of  the  PUFEM  spaces  is  the  notion  of  a  (M,  C^,  Cg )  partition  of  unity. 

Definition  1  Let  SI  C  Rn  be  an  open  set,  {fli}  be  an  open  cover  of  SI  satisfying  a 
pointwise  overlap  condition 

3 M  €  N  Wx  6  SI  card{i  |  x  6  fl;}  <  M. 
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Let  {c^}  be  a  Lipschitz  partition  of  unity  subordinate  to  the  cover  {ft,}  satisfying 


supp  ipi 

c 

closure(ft;)  Vi, 

(1) 

1  on  ft, 

(2) 

||¥5t||l,00(R") 

< 

C'oo, 

(3) 

||Vy>i||l,oo(Kn) 

< 

Cg 

(4) 

diam  ft;  ’ 

where  C Cg  are  two  constants.  Then  {tpt}  is  called  a  (M,  Coo,  Cg)  partition  of 
unity  subordinate  to  the  cover  {ft,} .  The  partition  of  unity  {95,}  is  said  to  be  of 
degree  ra  £  N0  if  {ipi}  C  Cm(R”).  The  covering  sets  {ft;}  are  called  patches. 

Definition  2  Let  {ft;}  be  an  open  cover  of  Q,  C  K"  and  let  {<£;}  be  a  (M,Coo,Cg) 
partition  of  unity  subordinate  to  {ft;}.  Let  Vi  C  l/1  (ft,  D  ft)  be  given.  Then  the  space 

V  ■=  ^  <PiVi  =  ipM  |  Vi  £  Vi}  C  H *(ft) 

i  i 

is  called  the  PUFEM  space.  The  PUFEM  space  V  is  said  to  be  of  degree  m  £  N  if 
V  C  The  spaces  V{  are  referred  to  as  the  local  approximation  spaces. 


Theorem  1  Let  Q  C  Rn  be  given.  Let  {tyi},  and  {Vi}  be  as  in  definitions  1, 

2.  Let  u  £  H 1(fi)  be  the  function  to  be  approximated.  Assume  that  the  local  approxi¬ 
mation  spaces  Vi  have  the  following  approximation  properties:  On  each  patch  0;  fl 
u  can  be  approximated  by  a  function  t \  £  Vi  such  that 

Ilu  - 'ui||i2(n,nn)  <  ei(i), 

||V(u  -  ,yi)||z,2(ninn)  <  e2(i). 


Then  the  function 

uap  =  Tivi  e  v  C  Hl{Sl) 

i 

satisfies 


n  ^ap||L2(n) 


||V(u-uap)||L2(n) 


< 


< 


vrMCoc 


) 


Proof:  Using  the  fact  that  (pi  =  1  on  3T2,  we  can  write  u  —  uap  =  J2i  <pi[u  —  vf). 
The  theorem  follows  after  an  application  of  the  second  estimate  of  lemma  2  with 

Ui  =  ipi(u  -  Vi).  □ 
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Example  1:  The  PUFEM  as  h  version.  Let  u  £  Hk(Q,),  k  >  1.  Let  each  patch 
have  diameter  hi  <  h,  and  let  each  Vi  have  approximation  properties 

ei(*)  <  C^i‘+1||wllH*(nnni))  /'c') 

e2(i)  < 

for  some  appropriate  p  >  0.  Then  the  error  estimates  of  theorem  1  take  the  form 

||u  -  ivIMn)  <  MCooW+^HIh^o) 

||V(«-tiap)||La(n)  <  MCy/2{CG  +  C00)h‘l\\u\\Hna)  {1 

where  we  used  the  first  estimate  of  lemma  2  in  the  estimate  of  the  sums  JT  ei(i)2, 
C2(z)2.  Note  that  estimate  (6)  holds  for  any  system  of  local  approximation  spaces 
Vi  satisfying  (5).  For  example,  if  the  spaces  Vi  consist  of  polynomials  of  degree  p,  then 
(5)  holds  with  p  =  min(A;  —  l,p).  If  the  spaces  Vi  consist  of  harmonic  polynomials  of 
degree  p,  (5)  holds  also  with  p  =  min(&  —  l,p)  if  we  know  a  priori  that  the  function  u 
is  harmonic.  In  this  example,  local  approximability  of  the  spaces  Vi  (and  thus  global 
approximability  by  theorem  1)  is  achieved  by  the  smallness  of  the  patches  0,  fl  Q. 

Example  2:  The  PUFEM  as  a  p  version.  Let  u  £  Hk(Q),  k  >  1,  and  let  {Di}£Li 
be  N  fixed  patches  covering  fi.  Denote  diam(f2i)  by  hi.  Assume  that  the  spaces  Vi 
(depending  on  a  parameter  p)  have  the  approximation  properties 

ei(z)  <  Chip-^Wullnk^nn.),  f?) 

e2(i)  <  C,p-M||u||Hfc(nnni)  1 

for  some  appropriate  p  >  0.  Then  the  error  estimates  of  theorem  1  take  the  form 


||«  -  Wop|U2(0)  <  M Coo C  max  hip  'i||«||jrk(n), 
l|V(u  -  U„„)|U,(n)  <  MC^2(ci  +  £?L)p-'‘||u||H.(fl|. 

Note  that  this  estimate  holds  for  any  system  Vi  satisfying  (7) — they  do  not  have  to 
polynomials  of  degree  p.  If  the  spaces  Vi  consist  of  polynomials  of  degree  p  then  (7) 
holds  with  p  =  k  —  1.  Estimate  (7)  also  holds  for  spaces  Vi  consisting  of  harmonic 
polynomials  of  degree  p  if  the  function  u  is  known  to  be  harmonic  (see  theorems  2, 
3).  In  this  example,  the  approximation  properties  of  the  global  PUFEM  space  are 
achieved  through  increased  approximability  of  the  local  spaces  while  keeping  the 
patches  fixed.  If  we  allow  the  size  of  the  patches  to  vary  as  well,  then  this  method 
behaves  like  an  hp  version. 


We  would  like  to  stress  at  this  point  that  the  requirements  on  the  partition  of  unity 
are  very  weak:  It  only  needs  to  be  Lipschitzian  in  order  to  produce  H 1  subspaces. 
Also,  we  do  not  need  positivity  of  the  partition  of  unity-the  elements  of  the  partition 
of  unity  are  allowed  to  change  sign.  However,  if  the  partition  of  unity  is  of  degree  m 
(and  the  local  approximation  spaces  are  sufficiently  smooth),  then  the  finite  element 


6 


space  V  as  constructed  in  definition  2  is  also  of  degree  m.  Theorem  1  is  formulated 
in  terms  of  H1,  appropriate  for  a  large  class  of  second  order  problems.  Mutatis 
mutandis  however,  the  estimates  can  be  formulated  in  terms  of  Hk,  k  >  1  to  produce 
finite  element  spaces  for  higher  order  equations.  Similar  estimates  can  be  achieved  in 
Sobolev  spaces  Wk,p. 

Remark  1:  This  idea  of  using  a  partition  of  unity  to  construct  finite  element 
spaces  tailored  to  the  differential  equation  has  been  used  in  [17],  [10],  and  [11],  As 
mentioned  in  the  introduction,  for  a  judicious  choice  of  parameters,  the  method  of 
[9]  reduces  to  a  special  type  of  PUFEM,  and  the  convergence  analysis  of  [13]  for  this 
special  case  is  based  on  theorem  1. 

3  The  PUFEM  in  One  Dimension 

3.1  A  One  Dimensional  Example 

Let  us  demonstrate  for  a  one  dimensional  model  problem  how  FE  spaces  with  good 
approximation  properties  are  constructed  with  the  PUFEM.  To  this  end,  consider 

-u"  +  k2u  =  feC2[ 0,1]  on  (0,1) 

u(0)  =  0  (8) 

■u'(l)  =  g  6  R. 

We  assume  that  the  parameter  k  >  1  is  large.  Associated  with  this  problem  is  an 
“energy”  norm,  given  by 

IMI-S  :=  {ll^lli^n)  +  P|Ml£2(n)} 

Let  us  note  that  for  large  k,  the  solution  to  problem  (8)  typically  exhibits  a  boundary 
layer  in  the  neighborhood  of  x  =  0,  and  thus  the  usual  FEM  perform  poorly  unless 
h  is  sufficiently  small  (relative  to  A:-1)  or  a  very  strongly  refined  mesh  is  used.  The 
PUFEM  allows  us  to  use  local  spaces  reflecting  this  behavior,  and  therefore  leads  us 
to  a  robust  FEM,  i.e.,  a  method  which  is  good  uniformly  in  k. 

Let  n  6  N,  h  =  L  and  define  Xj  =  jh,  j  =  0,...  ,n.  Define  also  x_!  =  —h, 
x„+i  =  1  +  h  and  let  the  patch  Qj  =  (xj_i, Zj+i),  j  =  0, ...  ,n.  On  each  patch  Q,j, 
we  have  to  define  a  local  space  which  can  approximate  the  solution  u  of  problem  (8) 
well.  We  consider 

Vj  =  span{l,  sinh  kx,  cosh  kx}  on  Q,j  fl  D,  j  =  l,...,n, 

Vq  =  span{sinh  kx,  1  —  cosh  kx}  on  Q0  0  Q,. 

We  note  that  the  space  Vo  is  constructed  such  that  it  satisfies  the  essential  boundary 
condition  at  x  =  0.  The  approximation  properties  of  these  spaces,  which  are  tailored 
to  this  particular  problem  (8),  are  given  by  the  following 
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Lemma  1  Let  u  be  the  solution  to  problem  (7)  and  let  flj,  V-  be  as  defined  above. 
Then  there  are  Vj  6  V-  such  that 


h2 , 


||(«-t*)#ll*>(n,nn)  <  Ch 1/2  /l2min(l,(^)-2)||/'||io0(n)  +  -||/licc(n) 
||(«  —  v»)llL*(n<nn)  <  Ch1/2  h3  min  (l,(kh)~2)  ||/'||i,oo(n)  + 


+y  A;  1)||/"||i,°o(n) 


where  C  >  0  is  independent  of  h,  k ,  and  f. 

Proof:  Because  the  spaces  V?  contain  the  fundamental  system  {sinh  kx,  cosh  kx },  it 
is  enough  to  approximate  a  particular  solution  to 

— u"  +  k2u  —  f  on  Q,j  fl  fl. 

By  Taylor’s  theorem,  on  fij  fl  fl,  f(x)  =  l(x)  +  r{x)  where  l{x)  is  linear  and  |r(x)|  < 
WII/l  i~(n)  (note  that  diam  Qj  <  2 h).  A  particular  solution  to  the  problem  with 
right  hand  side  r(x)  is  given  by  the  solution  ur  to 

—u"  +  k2ur  =  r  on  fl7-  fl  fl, 

u  =  0  on  5(flj  fl  fl). 


Thus, 


from  whence 


ltCllL2(njnn)  +  fc2||‘i*r||£a(ninn) 


2 

L°°(o,nn)> 


Kiu^no)  <  ch'/’f\\ru.ln), 

with  C  >  0  independent  of  h,  k,  and  /.  Finally,  a  particular  solution  to  the  problem 
with  right  hand  side  l(x)  is  given  by  ui(x)  =  k~2l(x )  which  can  be  approximated  in 
Vj  such  that 

\\ui  ~  vilU2(ninn)  +  h\\(ut  -  vJ)'||z,2(ninn)  <  C/i3/i1/2min(l,  (kh)~2)\\f\\L«,w, 

where  C  >  0  is  independent  of  h,  k,  and  /.  The  assertion  of  the  lemma  follows.  □ 

Remark  2:  The  spaces  Vj1  were  chosen  as  local  approximation  spaces  because 
they  contain  the  fundamental  system  {sinh  kx,  cosh  kx}  and  the  particular  solution 
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that  corresponds  to  a  constant  right  hand  side.  It  is  easy  to  check  that  the  func¬ 
tions  {1,  x, . . .  ,  xp}  actually  span  a  space  of  particular  solutions  for  polynomial  right 
hand  sides  of  degree  p.  Hence,  lemma  1  can  be  adapted  to  produce  the  following 
approximation  result.  The  spaces 

V?  =  span{sinh  kx,  cosh  kx,  1,  x, . . .  ,  xp}  on  Qj  fl  Cl,  j  =  l,...,n, 

Vq  =  span{sinh  kx,  1  —  cosh  kx,  x, . . .  ,  xp}  on  fi0  fl  Cl. 

contain  Vj  £  Vj  such  that 

||(u-t/i)'||1.(ninO)  <  Crh112  + 

+  ^ll/(f+2)lk-(0) 

||(u  —  Vi)||i2(ninn)  <  Cph1/2  [hp+3  min(l,{kh)  2)  ||/^p+1^||L“(n)+ 

up+ 2  1 

+  -j^-ram(h,  £)||/(p+2)||i,»(n) 

for  some  Cp  independent  of  h,  k ,  p,  and  /. 

For  any  partition  of  unity  {v?j}  subordinate  to  the  covering  {%},  the  finite  element 
space  V1  as  constructed  in  definition  2  is  given  by 

V1  =  span{y>j(z),  (pj(x)  sinh  kx,  tpj(x)  c°sh  kx, 

<Po{x)  sinh  kx,  </?0(®)(l  —  cosh  kx)  |  j  =  1, . . .  ,n}. 

Since  the  assumptions  on  the  partition  of  unity  stipulate  that  the  functions  (pj  be 
Lipschitz  continuous,  we  see  that  V1  C  Ff1(0).  Because  each  function  <pj  is  assumed 
to  vanish  outside  the  patch  Clj,  and  because  the  elements  of  V0r  vanish  at  x  =  0, 
we  see  that  all  elements  of  V1  vanish  at  x  =  0.  Hence,  a  conforming  finite  element 
method  can  be  based  on  V1 ,  and  the  finite  element  solution  is  the  best  approximant 
in  the  energy  norm: 

\\u-ufe\\e<  inf  IIu-vIIe. 

Therefore,  with  the  aid  of  theorem  1,  the  local  approximation  properties  of  the  spaces 
Vj  in  lemma  1  lead  to 

Proposition  1  Let  the  patches  {f2j}  and  the  local  approximation  spaces  {Vj1}  be 
given  as  above.  Let  {ipj}  be  a  (M,  ,  C’g)  partition  of  unity  subordinate  to  the 

patches  {fhj}.  Then  the  finite  element  solution  ufe  of  the  PUFEM  satisfies 

llu  —  ufe\\e  Ch2  {min  (l,  (kh)  1)  ||/,||Lco(n)  +  &  1||/,,||i,oo(n)}  (8) 


where  C  >  0  is  independent  of  h,  k,  and  f. 


This  shows  that  the  PUFEM  enables  us  to  construct  robust  finite  element  methods 
which  are  efficient  uniformly  in  k,  i.e.,  the  finite  element  method  behaves  as  well  for 
rough  case  of  large  k  as  it  does  for  the  smooth  case  k  —  1.  The  PUFEM  gives  these 
good  uniform  estimates  because  the  local  spaces  V?  capture  the  local  behavior  of  the 
exact  solution  very  well.  Note  that  the  number  of  degrees  of  freedom  is  comparable 
to  the  number  of  degrees  of  freedom  of  the  usual,  piecewise  quadratic  finite  element 
method  which  is  -  with  the  exception  of  (piecewise)  quadratic  solutions  -  of  order  h 2 
and  not  better.  Thus,  the  PUFEM  is  as  good  as  the  usual  piecewise  quadratic  finite 
element  method  for  the  smooth  case  k  =  1. 

A  simple  adaptation  of  this  idea  is  to  choose  the  local  spaces  selectively.  For  example, 
since  the  right  hand  side  /  is  smooth,  we  expect  a  boundary  layer  close  to  x  =  0  but 
expect  smooth  behavior  away  from  x  =  0.  Hence,  it  suffices  to  use  the  spaces  Vj  on 
patches  close  to  x  =  0,  and  we  can  use  polynomials  spaces  Vj  —  span{l,x, . . .  ,xp} 
on  patches  away  from  x  =  0.  The  idea  of  choosing  the  local  approximation  spaces 
selectively  can  also  be  employed  in  adaptive  versions  of  the  PUFEM.  Keeping  the 
patches  and  changing  the  degree  p  of  the  polynomials  lets  the  PUFEM  act  like  an 
adaptive  p  version;  changing  the  size  of  the  patches  adaptively  makes  the  PUFEM 
behave  like  an  adaptive  h  version. 

3.2  Examples  of  Partitions  of  Unity- 

In  this  section  we  propose  several  ( M ,  Coo,  Cg)  partitions  of  unity  for  the  one  dimen¬ 
sional  example  of  the  preceding  section.  Thus,  the  underlying  cover  of  the  domain 
(0, 1)  is  the  one  given  in  the  previous  section. 

1.  The  usual  piecewise  linear  hat-functions  form  a  partition  of  unity.  Let 

1  +  |  for  x  £  (— h,  0] 

1  -  |  for  x  G  (0,  h )  (9) 

0  elsewhere, 

and  define  the  partition  of  unity  by  <p){x)  =  <p(x  —  Xj),  j  =  0, . . .  ,n. 

2.  Functions  which  are  identically  1  on  a  subset  of  their  support  can  also  form  a 
partition  of  unity. 


<p\x)  =  l 


fi  +  2!  for  xe  (-|a,-j] 

2  J  1  for  x  €  (-£,  J) 

V(I)=  |-*t 

k  0  elsewhere, 

and  define  the  partition  of  unity  by  ifj(x)  —  (p(x  —  Xj),  j  =  0, . . .  ,n. 


(10) 
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3.  A  combination  of  the  above  two  examples  is  to  choose  the  functions  (p j  for 
patches  in  the  interior  but  to  modify  the  functions  on  patches  close  to  the 
boundary.  Define 

1  +  |  for  x  G  (— h,  0] 

1  for  x  G  (0,  h) 

2-  |  for  x  G  (h,  2h)  1  ’ 

0  elsewhere. 

We  observe  that  the  patches  D0  U  fix,  £ln-i  U  and  ttj,  j  =  2, . . .  ,  n  -  2,  cover 
D.  On  the  patches  f2  j,  j  =  2 , ...  ,n  —  2,  we  define  <p3(x)  =  y>)(x).  On 
patch  f2o  U  fii  we  choose  <p\(x)  =  <p3{x)  and  on  the  patch  fi„_i  U  we  choose 
v£-i(*)  =  fZ{.x  ~  Note  that  ip\  =  pi  +  ip\  and  p3n_x  =  <p\_x  + 

4.  In  all  three  examples  above,  the  partition  of  unity  is  merely  Lipschitz  continu¬ 
ous.  However,  partitions  of  unity  of  any  desired  regularity  can  be  constructed. 
Here  is  a  piecewise  polynomial  C 1  example.  The  resulting  global  finite  element 
space  Vn  is  then  a  subspace  of  C^O,  1].  Define 

(  (x  +  h)2(h  —  2x)  for  x  G  (—h,  0] 

<p4(x)  =  —  <  (h  —  x)2(h  +  2x )  for  x  G  (0,  h )  (12) 

( 0  elsewhere, 


and  define  the  individual  members  of  the  partition  of  unity  by  <p4(x)  =  (p4(x—Xj) 
on  the  patches  £lj. 


5. 


In  this  example,  let  Qj  be  any  cover  of  Q  satisfying  an  overlap  condition  (i.e., 
not  more  than  M  patches  overlap  in  any  given  point  igfl).  Let  ipj  be  Lipschitz 
continuous  functions  supported  by  the  patches  Qj.  If  <  C  and  »  — 

C'diamflj  on  each  flj  D  D,  for  some  C,  C  >  0  independent  of  j ,  then  the 
functions 


<Pj(X) 


(g) 

HiMx) 


form  a  (M,  CC-1,  (7C_1(1  +  MC2C~ 2))  partition  of  unity  subordinate  to  the 
cover  {flj}.  Note  that  the  functions  ipj  scale  with  their  supports  in  the  sense 
that  \ipj\  <  C  diamflj.  The  functions  cpj  inherit  the  smoothness  of  the  functions 
iftj,  i.e.,  with  this  “normalizing”  technique,  one  can  easily  construct  partitions 
of  unity  of  any  desired  regularity.  Another  feature  of  the  construction  is  that 
it  allows  us  to  build  (M,Coo,Cg)  partitions  of  unity  for  very  general  covering 
situations.  In  particular,  it  enables  us  to  produce  the  necessary  partitions  of 
unity  whenever  patches  are  added,  removed  or  otherwise  changed  in  an  adaptive 
computational  environment. 
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3.3  Polynomial  Local  Approximation  Spaces  and  Linear  De¬ 
pendencies 

In  this  section  we  want  to  analyze  in  more  detail  the  PUFEM  spaces  based  on  poly¬ 
nomial  local  approximation  spaces.  We  will  see  below  that  for  polynomial  local 
approximation  spaces,  the  choice  of  the  partition  of  unity  has  an  influence  on  the 
approximation  properties  of  the  PUFEM  space  and  has  implement  ational  ramifica¬ 
tions  in  the  following  sense.  In  any  implementation,  a  basis  of  the  PUFEM  space 
has  to  be  constructed  and  it  would  convenient  if  that  basis  were  determined  directly 
by  the  basis  functions  of  the  local  approximation  spaces.  In  general  however,  this  is 
not  true.  For  example,  for  piecewise  linear  partitions  of  unity  and  polynomial  local 
approximation  spaces,  the  local  basis  functions  (multiplied  by  appropriate  partition 
of  unity  function)  are  linearly  dependent  and  thus  do  not  form  a  basis  of  the  PUFEM 
space  (see  below).  Although  this  example  is  artificial,  it  suggests  that  even  if  the  local 
basis  functions  lead  to  a  basis  of  the  PUFEM  space,  the  resulting  functions  might  be 
“nearly”  linearly  dependent,  and  the  resulting  finite  element  stiffness  matrix  will  be 
badly  conditioned. 

Define 

V0n  =  span{x, . . .  ,  xp}  on  fl0  H  fl, 

Vj1  =  span{l,  x, . . .  ,  xp}  =  span{l,  x  —  Xj, . . .  ,  (x  —  Xj)p }  on  f lj  fl  fl 

for  j  =  1, . . .  ,n  and  set  Vj1  =  {0}  if  p  =  0.  For  any  partition  of  unity  {<Pj},  the 
PUFEM  space  is  given  by 

V11  =  spa.n{(pj(x)xm,<p0(x)xq\j  =  1,...  ,n,  m  =  0,...,p,  q  =  l,...,p}. 

The  fact  that  the  functions  {<pj}  form  a  partition  of  unity,  i.e.,  Ylj(Pj(x)  =  1  on 
fl,  implies  that  the  space  V11  satisfies  a  consistency  condition  in  the  sense  that  all 
polynomials  of  degree  <  p  which  vanish  in  x  =  0  are  contained  in  Vn . 

Let  us  now  consider  the  spaces  V11  based  on  the  various  partitions  of  unity  of  the 
previous  section  more  closely.  Denote  by  V11’1,  V11’2,  and  V11,3  the  spaces  Vn  con¬ 
structed  using  the  partitions  of  unity  {y^},  and  {v’j}  respectively.  Let  us 

concentrate  on  V11’1  first.  Owing  to  the  fact  that  the  functions  cp)  are  piecewise 
polynomials,  the  space  V11,1  is  precisely  the  space  of  piecewise  polynomials  of  degree 
p  +  1  constrained  to  vanish  in  x  =  0.  This  is  an  example  where  the  global  finite  ele¬ 
ment  space  has  even  better  approximation  properties  than  guaranteed  by  theorem  1: 
Locally,  approximation  is  done  by  polynomials  of  degree  p  and  theorem  1  states  that 
the  local  approximation  properties  are  inherited  by  the  global  space,  i.e.,  the  H 1  ap- 
proximability  is  0(hp).  However,  the  space  of  piecewise  polynomials  of  degree  p  +  1 
has  better  approximation  properties:  it  is  0(hp+1)  for  H1  estimates.  Let  us  note  that 
dimU^7,1  =  n(p  +1). 
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As  mentioned  above,  it  would  be  convenient  for  implementational  purposes  to  take  as 
a  basis  of  the  finite  element  space  V11'1  functions  which  are  determined  by  the  basis 
functions  of  the  local  spaces  V?1,  i.e.,  we  would  like  to  take  the  functions 

<p)(x)(x  -  Xj)m  y  j  =  1, ...  ,n,  m  =  0,...,p,  (13) 

<pj(x)xm,  m  =  l,...,p.  (14) 

However,  these  functions  are  not  linearly  independent  for  p  >  1  as  a  simple  counting 
argument  reveals:  there  are  n(p+l)+p  functions  but  dim  V11,1  =  n(p+l)  <  n(p-f-l)+p 
for  p  >  1.  Of  course,  one  can  still  use  these  functions.  For  problem  (7)  they  will  lead 
to  a  positive  semi-definite  matrix  (as  opposed  to  a  positive  definite  matrix,  which  is 
obtained  if  a  basis  is  used),  which  has  many  algebraic  solutions.  However,  all  these 
algebraic  solutions  are  merely  different  representations  of  the  same  function  on  Q. 
One  way  to  solve  this  linear  system  is  to  use  a  penalty  method  to  deal  with  the  linear 
dependencies  (see  [10]  for  a  computational  analysis). 

One  can  avoid  these  linear  dependencies  if  one  uses  a  different  partition  of  unity.  For 
example,  whenever  the  partition  of  unity  is  such  that  each  member  (pj  is  identically  1 
on  an  open  set  Oj  C  flj  fl  Q  (and  all  the  other  ones  vanish  there),  linear  dependencies 
as  above  cannot  occur.  Hence,  the  functions 

<p!(x)(x  -  Xj)m,  j  =  l,...,n,  m  =  0,...,p,  (15) 

cpo(x)xm,  m  =  l,...,p.  (16) 

form  indeed  a  basis  of  the  space  V11,2. 

A  more  careful  analysis  of  the  linear  dependencies  occurring  for  the  case  of  V11’1 
reveals  that  the  local  approximation  space  at  either  the  left  or  the  right  endpoint 
of  fi  contains  too  many  functions.  Thus,  a  modification  of  the  partition  of  unity 
at  one  (or  both)  endpoints  allows  us  to  exclude  linear  dependencies:  The  functions 
<Pj(x)(x  —  Xj)m,  j  =  2, . . .  ,n  —  1,  m  =  0, . . .  ,p,  <p3(x)xm,  m  =  1, . . .  ,p,  form  a  basis 
of  V11,3.  Let  us  point  out  that  the  space  V11,3  does  no  longer  contain  all  piecewise 
polynomials  of  degree  p  +  1.  Let  us  note  here  that  this  space  is  very  closely  related 
to  V11,1.  In  fact  for  problem  (7)  the  stiffness  matrix  of  the  finite  element  method 
based  on  V11  ,3  can  be  easily  extracted  from  the  positive  semi-definite  stiffness  matrix 
constructed  using  V11’1. 

The  example  V11,1  shows  that  “unfortunate”  combinations  of  local  approximation 
spaces  and  partitions  of  unity  exist,  where  the  basis  elements  of  the  local  spaces 
multiplied  by  the  appropriate  partition  of  unity  function  are  linearly  dependent.  This 
indicates  that  even  if  the  chosen  functions  derived  from  the  local  bases  are  linearly 
independent  and  form  a  basis  of  the  finite  element  space,  the  resulting  stiffness  matrix 
may  still  be  badly  conditioned. 
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3.4  Polynomial  Local  Approximation  Spaces;  Lagrange  Type 
Elements 

If  we  choose  the  functions  <pj(x)(x  —  Xj)m  as  basis  functions  of  the  space  V11,  the 
degrees  of  freedom  cannot  be  identified  directly  as  function  values  in  certain  points. 
Rather,  the  degrees  of  freedom  are  related  to  higher  derivatives  of  the  elements  of 
Vn  in  the  points  Xj.  In  this  sense,  the  functions  ipj(x)(x  —  Xj)m  produce  a  Hermite 
type  space.  However,  it  is  also  possible  to  construct  Lagrange  type  spaces,  where  the 
degrees  of  freedom  represent  the  function  values  in  particular  “Lagrange  interpolation 
points”.  Let  us  illustrate  this  for  the  case  where  we  want  to  approximate  locally 
with  polynomials  of  degree  p.  Let  {flj}  be  a  cover  of  £2  =  (0, 1)  and  let  {<pj}  be  a 
(M,  Coo,  Cq)  partition  of  unity  subordinate  to  the  cover.  Let  yi,  i  =  1, . . .  ,  N,  be  the 
“Lagrange  interpolation  points” ,  and  assume  that  there  are  p  +  1  points  yi  in  each 
patch  flj.  In  order  to  be  able  to  enforce  the  essential  boundary  condition  at  x  =  0,  we 
will  stipulate  y\  =  0.  On  each  patch  fij,  let  LjtVi  be  the  usual  polynomial  Lagrange 
interpolation  function  of  degree  p  which  is  1  in  the  point  y,  and  vanishes  in  all  the 
other  p  “Lagrange  interpolation  points”  which  are  in  the  patch  Qj.  As  before,  we 
define  the  global  finite  element  space  by 


~  I  ^ j  ( x  )  x )  aj,yi  I  aj,y,  £ 


This  is  exactly  the  same  space  as  is  obtained  if  the  local  spaces  Vj  are  chosen  to  be 
span{l,x, . . .  ,  xp}.  Now,  if  we  identify  unknowns  associated  with  the  same  interpo¬ 
lation  point,  i.e.,  if  we  set  am>v;  =  a„)Vi  for  all  n,  m  for  each  point  yi,  and  denote  these 
common  values  by  ayi,  we  arrive  at  the  space 


V  —  <  'y  '  <Pj{X)Lj,yi(x)  ay;  ayx  £  R.  ’  • 

[  t=l  -I 

Because  the  functions  <pj  form  a  partition  of  unity  and  because  the  functions  LjiVi 
take  only  the  values  0  and  1  in  the  “Lagrange  interpolation  points”  ym,  the  values 
ayi  are  precisely  the  function  values  of  the  elements  of  VIv.  Hence,  we  can  take  as  a 
basis  of  Viv  the  functions 


$,(x)  —  (pj(x)Lj>yi(x), 


i  =  !,•••  ,N. 


The  essential  boundary  condition  at  x  =  0  is  also  easily  enforced  by  simply  setting 
ayi  =  0,  which  gives  the  space 


14 


Let  us  make  a  few  remarks  on  the  approximation  properties  of  the  space  Vv .  The 
approximation  properties  of  the  spaces  V 11  are  given  by  the  approximation  proper¬ 
ties  of  the  local  spaces  V-11,  i.e.,  by  the  approximation  properties  of  polynomials  of 
degree  p.  For  fixed  degree  p  and  appropriate  conditions  on  the  distributions  of  the 
interpolation  points  on  each  patch,  it  can  be  shown  that  approximation  with  Vv  is 

-  up  to  a  constant  -  as  good  as  with  V11 .  Finally,  let  us  mention  that  these  spaces 

y/V,  yV 

are  closely  related  to  the  method  proposed  in  [9]. 

Above  we  noted  that  the  close  relation  between  the  spaces  V11,1  and  V11,3  enables  us 
to  construct  the  stiffness  matrix  based  on  V11,3  easily  from  the  one  based  on  V11’1. 
Similarly,  the  stiffness  matrix  based  on  the  functions  <f>,  can  be  extracted  from  the 
stiffness  matrix  based  on  the  functions  <pjLjiy{. 

4  Comments  on  Choosing  Local  Approximation 
Spaces 

4.1  Change  of  Variables  Techniques 

In  section  3.1,  we  chose  the  local  approximation  spaces  for  problem  (7)  to  consist 
of  a  fundamental  system  for  the  differential  equation  and  particular  solutions  for 
polynomial  right  hand  sides.  A  different  method  to  construct  local  approximation 
spaces  is  based  on  changes  of  variables.  If  the  change  of  variables  x  t— >  x  maps  the 
problem  onto  a  problem  which  can  be  approximated  well  (in  some  appropriate  norm) 
by  polynomials  (in  x),  say,  then  the  “mapped  polynomials”,  i.e.,  P(x(x)),  where  P  is 
a  polynomial,  also  have  good  approximation  properties.  For  example,  [17]  considered 
the  problem 

—dx(a(x,y)dxu)  —  dy(a(x,y)dyu)  =  /  on  Q, 

u  =  0  on  dCl 

where  the  coefficient  a{x,y)  is  assumed  to  satisfy 

0  <  a  <  a(x,  y)  <  /3  <  oo 

and  is  uni-directionally  rough,  i.e.,  the  coefficient  a(x,y )  is  smooth  in  the  y  direction 
while  it  is  rough  in  the  x  direction.  The  roughness  of  the  coefficient  a(x,y )  results  in 
poor  regularity  properties  of  the  solution  u,  and  thus  the  usual  finite  element  method 
leads  to  mesh  sizes  h  which  are  prohibitively  expensive.  For  the  simplified  model, 
a(x,y )  =  a(x),  the  change  of  variables 


transforms  the  problem  into  one  for  which  a  better  regularity  theorem  holds:  if  f  £ 
L2(£l),  then  the  transformed  function  u  is  in  H2(£l)  (fi  denotes  the  image  of  Q,  under 
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the  above  transformation;  cf.  [17]  for  a  proof).  Thus,  u  can  be  approximated  by 
linear  functions  in  x,  y.  Formulating  in  the  original  coordinates  gives  that  u  can  be 
approximated  on  the  patch  Qj  by 


Vj  € 


dt 

Wr 


y } 


such  that 


11“  “  ^  C  (diamfij)  ||/ 

The  constant  C  >  0  depends  only  on  a,  (3  and  is  independent  of  the  roughness  of  the 
coefficient  a(x)  and  thus  these  local  spaces  have  good  approximation  properties  on 
patches  independent  of  the  bad  behavior  the  coefficient  a(x)  might  display. 

Let  us  finally  point  out  for  this  example  that  the  change  of  variables  can  be  done 
locally:  if  (xj,yj)  £  flJ}  then  the  approximating  functions  can  be  chosen  to  be  in 
sPan{l,/^  y}. 

Another  instance  where  the  idea  of  using  a  change  of  variables  is  successfully  used 
can  be  found  in  [14,  15].  For  elliptic  problems  in  2-d  with  corners  or  interfaces,  the 
use  of  a  conformal  map  is  proposed  which  maps  the  rough  solution  to  a  smoother 
function.  This  smoother  function  on  the  mapped  domain  can  be  approximated  by 
polynomials.  Hence,  the  images  of  polynomials  under  the  inverse  of  this  conformal 
map  are  used  for  the  approximation  of  the  original  problem. 


4.2  Optimality  of  Local  Approximation  Spaces  and  n— Width 

An  interesting  issue  in  the  context  of  finding  good  local  approximation  spaces  is  the 
question  of  optimality  of  the  local  spaces.  We  measure  optimality  in  terms  of  n- width, 
i.e.,  in  terms  error  per  degree  of  freedom  for  a  whole  class  of  functions: 

d(n,  II  •  II.  S)  =  inf  sup  inf  ||/  -  g\\ 

xin  f£S  9^-^tl 

where  En  denotes  an  n-dimensional  space,  and  S  is  the  class  of  functions  that  we  wish 
to  approximate;  typically,  S  is  chosen  as  the  unit  ball  of  some  appropriate  Banach 
space.  A  minimizing  space  En  is  called  an  optimal  space.  We  see  that  this  notion 
of  optimality  depends  on  the  dimension  n,  the  norm  ||  •  ||,  in  which  we  measure  the 
approximation  error,  and  the  choice  of  the  class  S.  In  particular,  different  classes  S 
lead  to  different  optimal  spaces.  In  practice,  of  course,  we  want  robust  optimal  (or 
near  optimal)  approximation  spaces  because  we  might  not  know  with  respect  to  which 
class  of  functions  we  should  optimize  (this  uncertainty  issue  is  elaborated  in  [16]).  For 
example,  if  we  choose  ||  •  ||  =  ||  •  jjiyi(n) j  and  if  we  are  interested  in  approximating 
functions  which  are  analytic  on  Cl  DD  fi,  the  class  S  could  be  taken  as  the  unit  ball  of 
any  Hk(Q),  k  >  1.  Thus,  due  to  this  uncertainty,  we  want  the  approximation  spaces 
to  be  optimal  for  as  large  a  class  of  functions  as  possible. 
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Proposition  2  below  exhibits  an  example  of  approximation  spaces  which  are  optimal 
for  large  classes  of  harmonic  functions.  In  the  framework  of  the  PUFEM,  proposi¬ 
tion  2  yields  the  following  result:  For  the  approximation  of  harmonic  functions  on 
patches  fi,  which  are  discs  C  M2,  the  choice  of  spaces  of  harmonic  polynomials  as 
local  approximation  spaces  V{  is  the  optimal  one.  The  notion  of  optimality  here  is 
tied  to  the  assumptions  of  proposition  2,  namely,  the  restriction  to  functions  defined 
on  discs  and  to  rotationally  invariant  classes  of  harmonic  functions. 

For  ease  of  exposition,  we  deal  with  complex-valued,  holomorphic  (analytic)  functions 
and  observe  that  the  case  of  harmonic  functions  follows  by  taking  real  parts. 

We  introduce  the  spaces 

Hk  =  {/  €  Hk{Bn{ 0))  |  /  is  holomorphic  on  .Br(O)},  k  >  0 

and  a  Hilbert  space  Ti  of  holomorphic  functions  with  inner  product  <  -,  •  >r  and 
norm  ||  •  \\ji  is  called  rotationally  invariant  if 

ll/Mllw  =  Il/(*)||w  Vae<C,M  =  l. 

The  space  7i°  (771)  is  a  Hilbert  space  with  the  L2  (H1)  inner  product  and  thus  a  closed 
subspaces  of  L2(Br( 0))  (H'^Br^O))).  Therefore,  the  space  L2(Br( 0))  (B1(Br(0))) 
can  be  written  as  the  direct  sum  of  'H0  (H1)  and  its  orthogonal  complement.  This 
reduces  the  search  for  optimal  spaces  for  the  approximation  of  holomorphic  functions 
in  the  L2  (H1)  norm  to  the  problem  of  finding  optimal  subspaces  of  TiP  (H1). 

The  polynomials  (.zn)£L0  form  an  orthogonal  basis  of  'Hk,  k  >  0,  and  it  is  easy  to 
see  that  they  actually  form  an  orthogonal  basis  for  any  rotationally  invariant  Hilbert 
space  of  holomorphic  functions.  Therefore,  setting  V’«(n)  =<  zn,zn  >r,  gives  the 
representation 


ll/llw  =  i/"i2^(n) 

n=0 

where  the  fn  are  the  Taylor  coefficients  of  the  holomorphic  function  /,  i.e., 


f(z)  =  -f"2"  on 

71—0 


For  example,  we  have 

if)no(n)  =  7 T—^—R2n+2, 
n  +  1 

V>wi(n)  =  TrR2n  +  nj  . 
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Proposition  2  Let  Hi,  H2  be  two  rotationally  invariant  Hilbert  spaces  of  holomor- 
phic  functions  on  Br( 0).  Assume  that  the  quotient 


TpH2(n) 


is  monotonically  decreasing  in  n.  Then  the  spaces 


Tp  =  span{zn  |  n  =  0, . . .  ,p} 


are  optimal  spaces  for  the  approximation  of  functions  in  H2  in  the  ||  •  \\'ul  norm,  i.e., 
the  space  Tp  minimizes  the  expression 


sup  inf 
fen2  9eEr 


11/  ~  9b<i 

\\f\\n2 


(17) 


over  all  p  dimensional  subspaces  Ev  of  Hi . 

Proof:  The  proof  proceeds  in  two  steps.  First,  we  will  see  that  (17)  is  bigger  than 
or  equal  to 

/  ^Mp  +  i)\1/2 

Ww2(p  +  i)y 

for  any  p  dimensional  subspace  of  T~Ci.  In  the  second  step,  we  see  that  this  infimum 
is  attained  for  the  choice  of  Tp  as  p  dimensional  approximation  space.  Let  a  p  di¬ 
mensional  subspace  Ev  of  Tix  be  given.  Choose  /  £  Tp+i  orthogonal  (with  respect  to 
<  •)  •  >7ii)  to  Ep.  Then,  the  square  of  (17)  can  be  bounded  from  below  by 

ll/llt,  '  !*1  =  Wp+i] 

ll/llt,  -  ll/ll?,.  M  +  i) 

where  we  made  use  of  the  monotonicity  assumption.  On  the  other  hand,  the  choice 
Ep  =  Tp  implies 

•  {  \\f-9Wn,  ^  Er=P+i  l/nl ^  ipni(p  +  1) 

||/||^  -JZ  E“  „  l/.IVx.W  -^,(P+1) 

where  we  made  again  use  of  the  monotonicity  assumption.  □ 

Choosing  TL\  in  proposition  2  to  be  H°  or  Hl  shows  that  the  spaces  Tp  axe  optimal  if 
we  measure  approximability  in  the  L2  or  H1  norm  and  if  we  approximate  rotationally 
invariant  classes  which  satisfy  a  certain  monotonicity  of  the  numbers  ip-H2 (n).  All 
spaces  Hk  fall  into  this  latter  category  and  many  spaces  of  holomorphic  functions 
which  are  in  some  weighted  Hk  spaces  where  the  weight  is  rotationally  symmetric. 
Let  us  further  note  that  in  the  context  of  the  PUFEM,  theorem  1  suggests  that  we 
optimize  with  respect  to  the  norm  (diam2(fii)|  •  +  II  '  llL’(no)1/2-  The  Proof  of 
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proposition  2  shows  that  this  choice  of  norm  also  leads  to  the  spaces  Tp  as  optimal 
approximation  spaces. 

Remark  3:  As  stated  earlier,  proposition  2  can  be  formulated  for  harmonic  func¬ 
tions  as  well.  Then,  the  2p  +  1  dimensional  spaces  of  harmonic  polynomials  are 
optimal  under  similar  conditions.  For  example,  the  2p  +  1  dimensional  spaces  of  har¬ 
monic  polynomials  are  optimal  for  the  approximation  of  harmonic  function  on  the 
discs  jBr(O)  which  are  in  the  spaces  Hk(BR( 0)),  k  >  1. 

Remark  4:  Proposition  2  and  the  preceding  remark  state  (loosely  speaking)  that 
harmonic  polynomials  are  universally  optimal  for  the  approximation  of  harmonic 
functions  on  discs.  This  is  partly  a  justification  for  the  approximation  with  harmonic 
polynomials  in  section  6.1:  As  long  as  the  patches  differ  not  too  much  from  discs,  we 
expect  spaces  of  harmonic  polynomials  to  be  nearly  optimal  for  the  approximation  of 
harmonic  functions. 

Let  us  stress  here  that  harmonic  polynomials  are  no  longer  optimal  if  one  of  the 
assumptions  of  proposition  2  is  changed.  For  example,  consider  approximation  on  a 
sector  W  with  angle  u>  and  size  R  (for  notational  convenience,  we  identify  R2  with 
the  complex  plane  C): 

W  —  {z  E  C  |  \z\  <  R  and  0  <  axg  z  <  w}. 

Assume  that  we  are  interested  in  approximating  (in  H1,  say)  harmonic  functions 
satisfying  homogeneous  Dirichlet  conditions  on  the  two  straight  sides  of  the  sector, 
i.e.,  functions  of  the  form 


OO 

u  =  ^2anImzm/u 

Tl—l 


with  coefficients  an  E  R.  Then  the  functions  Imz™/",  n  =  1, ...  ,p  form  optimal 
spaces  of  dimension  p  for  the  whole  scale  of  spaces 


7? 


^anImzm/“ 

71=1 


anE  R  and  |an|2(l  +  nf^R2^  < 

n= 1 


k  >  1. 


The  proof  of  this  statement  is  very  similar  to  the  proof  of  proposition  2.  A  different 
way  of  defining  the  spaces  fft  is  to  say  that  harmonic  functions  in 
which  are  antisymmetric  with  respect  to  the  real  axis  are  mapped  onto  the  elements 
of  Hk  under  the  conformal  change  of  variables  z  i-»  W*. 


5  The  PUFEM  in  Two  Dimensions 

In  the  two  dimensional  case  -  just  as  in  the  one  dimensional  one  -  we  have  to  address 
the  creation  of  a  partition  of  unity  and  the  choice  of  local  approximation  spaces.  Let 
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us  first  outline  two  different  types  of  partitions  of  unity.  If  a  domain  D  fl  is  given 
by  a  mesh  (i.e.,  triangles,  quadrilaterals,  or  mapped  triangles  or  quadrilaterals),  then 
the  usual  pyramid  functions  associated  with  the  nodes  of  the  mesh  form  a  piecewise 
smooth  partition  of  unity.  Since  in  all  the  numerical  examples  below,  we  use  this  kind 
of  partition  of  unity,  let  us  exemplify  this  idea  with  one  example.  Let  Q  be  the  unit 
square  (0, 1)  x  (0, 1)  and  let  it  be  subdivided  into  n2,  n  6  N  squares  of  side  length 
h  =  -  with  nodes  ( Xj,yj ),  j  =  1, . . .  ,  (n  +  l)2.  Define 


<p(x) 


'{1  -x)(l -y) 

(1  +x)(l  -  y) 

<  (1  +  x)(l  +y ) 

(1  -x)(l  +y) 

0  elsewhere. 


for  (x,y)  G  [0,1]  x  [0,1] 
for  ( x,y )  6  [-1,0]  x  [0,1] 
for  (x,  y )  e  [-1,0]  x  [-1,0] 
for  (x,y)  e  [0,1]  x  [-1,0] 


(18) 


Then  the  functions  tpj(x)  =  y?((x  —  Xj)/h,(y  —  yj)/h )  associated  with  the  (n  +  l)2 
patches  Ctj  =  {(x,y)  |  |x  —  Xj\  <  h,  \y  —  yj\  <  h}  form  a  partition  of  unity.  This  is  the 
analogous  construction  to  the  first  construction  of  section  3.2. 

The  second  type  of  partition  of  unity  is  given  by  the  construction  described  in  the 
fifth  method  of  section  3.2.  For  example,  if  fi  is  covered  by  circles,  ellipses,  or  quadri¬ 
laterals,  it  is  easy  to  construct  a  partition  of  unity  of  any  desired  regularity  by  the 
“normalizing”  technique  outlined  in  the  fifth  method  of  section  3.2. 

Let  us  stress  at  this  point  that  the  partition  of  unity  does  not  have  to  be  related  to 
the  geometry  of  the  domain  of  interest. 

Many  of  the  observations  of  section  3.3  about  the  one  dimensional  case  are  true  in 
the  two  dimensional  setting  as  well.  For  example,  it  can  be  shown  that  the  piecewise 
bilinear  partition  of  unity  described  above  in  conjunction  with  polynomial  local  ap¬ 
proximation  spaces  Vj  displays  the  same  difficulties  with  linear  dependencies  as  the 
space  V11,1  of  section  3.3  (cf.  [10]).  However,  the  same  idea  of  modifying  the  parti¬ 
tion  of  unity  on  patches  close  to  the  boundary  as  is  proposed  in  the  third  method  of 
section  3.2  leads  to  a  basis  of  the  finite  element  space  which  is  directly  related  to  the 
bases  of  the  local  spaces.  As  observed  in  the  one  dimensional  case,  the  stiffness  matrix 
resulting  from  the  modified  partition  of  unity  can  actually  be  constructed  from  the 
original  one. 

Related  to  the  choice  of  the  partition  of  unity  (and  of  the  local  approximation  spaces) 
is  the  question  of  integrating  the  shape  functions  against  each  other,  because  the  par¬ 
tition  of  unity  is  typically  only  piecewise  smooth  (and  hence  the  shape  functions). 
This  issue  will  be  explored  in  more  details  in  a  forthcoming  paper.  For  all  the  nu¬ 
merical  examples  below,  we  use  the  partition  of  unity  for  the  unit  square  described 
above,  and  therefore  the  usual  integration  schemes  on  each  of  the  n2  square  can  be 
applied. 

Another  important  question  is  the  implementation  of  essential  boundary  conditions. 
For  some  problems,  it  is  easy  to  construct  local  approximation  spaces  Vj  on  patches 
close  to  the  boundary  which  have  both  good  approximation  properties  and  satisfy  the 
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essential  boundary  conditions.  This  is  the  case  in  the  one  dimensional  problem  (7) 
with  the  choice  Vj .  For  an  example  in  two  dimension,  consider  the  implementation  of 
homogeneous  Dirichlet  conditions  on  a  straight  part  of  the  boundary  for  the  problem 
—  A u  =  0.  Here,  harmonic  polynomials  which  are  antisymmetric  with  respect  to 
that  straight  line  have  good  approximation  properties  and  satisfy  the  homogeneous 
boundary  conditions.  A  similar  approach  works  in  a  corner. 

One  way  to  imitate  the  way  essential  boundary  conditions  are  implemented  in  the 
classical  finite  element  methods  is  to  use  spaces  of  (piecewise)  full  polynomials  on 
patches  close  to  the  boundary.  In  that  case,  all  the  techniques  of  the  usual  finite 
element  methods  can  be  applied.  Another  approach  to  the  implementation  of  essential 
boundary  conditions  is  the  use  of  Lagrange  multipliers  or  a  penalty  method.  In  the 
numerical  examples  below,  we  chose  the  boundary  conditions  to  be  natural  in  order 
to  be  able  to  concentrate  on  the  approximation  properties  of  the  spaces  constructed 
with  the  PUFEM. 

6  Numerical  Examples 

In  this  section,  we  will  present  two  numerical  examples,  namely,  the  approximation 
of  solutions  to  Laplace’s  equation  and  Helmholtz’s  equation  on  the  unit  square  with 
the  PUFEM. 

6.1  Laplace’s  Equation 

Let  us  consider  first  approximations  to  the  solution  of 
—A u  =  0  on  V,  =  (0, 1)  x  (0, 1) 
dnu  =  <9nRe(-J— ^  +  ^1—)  on  dQ,  a  =1.05, 

and  we  fix  u  in  (0, 0)  in  order  to  make  the  solution  of  this  problem  unique.  Since 
we  want  to  present  a  p  version  of  the  PUFEM  where  the  local  approximation  spaces 
are  chosen  as  spaces  of  harmonic  polynomials  of  degree  p,  we  need  to  clarify  the 
approximation  properties  of  harmonic  polynomials.  This  is  done  in  the  following  two 
theorems.  Note  that  there  are  only  2p  +  1  harmonic  polynomials  of  degree  p. 

Theorem  2  (Szego)  Let  VI  C  R2  be  a  simply  connected ,  bounded  Lipschitz  domain. 
Let  fi  DD  Q  and  assume  that  u  £  L2(£l)  is  harmonic  on  (7.  Then  there  is  a  sequence 
(up)£L0  of  harmonic  polynomials  of  degree  p  such  that 

||u  -  UplU-(n)  <  C'e-7I’||u||i2(Q), 
l|V(«-Up)|U-(n)  <  Ce  7p|M|La(fi) 

where  7,  C  >  0  depend  only  on  VI,  VI. 


Proof:  See  [21],  [24]. 


□ 


Theorem  3  Let  Cl  be  a  simply  connected  bounded  Lipschitz  domain,  star-shaped  with 
respect  to  a  ball.  Let  the  exterior  angle  of  Cl  be  bounded  from  below  by  Xtt,  0  <  A  <  2. 
Assume  that  u  €  Hk{Cl),  k  >  1,  is  harmonic.  Then  there  is  a  sequence  (up)£L2  °f 
harmonic  polynomials  of  degree  p  such  that 

1  \  X(k :~j) 

j  j  = 

where  C  >  0  depends  only  on  Cl  and  k. 

See  [11]  for  a  proof  of  theorem  3.  Note  that  typically  A  <  1  and  that  for  domains 
with  re-entrant  corners,  A  can  be  significantly  less  than  1. 

Remark  5:  The  restriction  in  theorem  3  that  Cl  be  star-shaped  with  respect  to 
a  ball  is  not  a  big  constraint  for  our  purposes  because  we  are  interested  in  local 
estimates  on  patches  and  the  patches  are  typically  chosen  to  be  star-shaped. 


u  -  Up||fri(n)  <  C 


Figure  1:  PUFEM,  classical  p  version  for  Laplace’s  equation;  a  =  1.05, 8x8  elements 


For  the  PUFEM  the  domain  Cl  is  covered  by  square  patches  and  the  partition  of  unity 
is  chosen  to  be  piecewise  bilinear  as  described  in  section  5.  The  specific  choice  n  =  8 
was  made,  and  the  local  approximation  spaces  Vj  consist  of  harmonic  polynomials 
of  degree  p  (p  ranging  from  0  to  8).  In  figure  1  we  plot  the  relative  error  in  energy 
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norm  (i.e.,  the  relative  error  in  the  H1  semi  norm)  versus  the  number  of  unknowns 
for  three  methods.  The  PUFEM  is  compared  with  two  classical  p  versions,  namely, 
the  tensor  product  spaces  Qp  and  the  serendipity  spaces  Q'  based  on  an  8  x  8  mesh. 
We  see  clearly  that  the  use  of  harmonic  polynomials  made  possible  by  the  PUFEM 
is  advantageous:  in  order  to  achieve  1%  error  in  H1,  the  PUFEM  based  on  harmonic 
polynomials  needs  only  half  as  many  DOF  as  the  usual  p  version  spaces  Qp,  Q'p.  This  is 
in  accordance  with  our  earlier  observation  that  the  number  of  harmonic  polynomials 
grows  linearly  with  the  degree  p  whereas  the  size  of  full  polynomial  spaces  grows 
quadratically.  Note  that  the  disparity  between  the  PUFEM  and  the  spaces  of  full 
polynomials  becomes  bigger  for  higher  accuracy.  See  [10]  for  a  more  detailed  study 
of  the  performance  of  the  PUFEM  as  the  parameters  n  and  a  are  varied. 

Remark  6:  For  the  elasticity  equations  in  2-d,  the  situation  is  completely  analogous 
to  Laplace’s  equation.  In  the  absence  of  body  forces,  the  displacement  field  (u,v) 
under  the  plane  strain  assumption  can  be  expressed  by  two  holomorphic  functions  ip , 
ip  (see  [12]): 

2p(u  +  iv)  =  K(p(z)  —  z<p'(z)  —  ip{z)  (19) 

where  k  =  (A  +  3/x) / (A  +  p)  and  A,  p  are  the  Lame  constants.  Choosing  k  = 

(A*  +  3p)/(\*  +  p)  with  A*  =  2A/i/(A  +  2 p)  gives  the  representation  for  the  case 

of  plane  stress.  The  holomorphic  functions  (p ,  ip  can  be  approximated  by  complex 
polynomials  tpp,  ipp  of  degree  p,  and  thus  the  functions 

K<pp(z)  -  zlpp(z )  -  ipp(z) 

take  the  role  of  “harmonic”  polynomials  for  the  elasticity  equations.  It  can  be  shown 
that  theorems  2,  3  hold  verbatim  for  the  approximation  of  the  solutions  to  the  elas¬ 
ticity  equations  with  these  “harmonic”  polynomials  (see  [11]). 

6.2  Helmholtz’s  Equation 

The  next  numerical  example  deals  with  the  approximations  to  Helmholtz’s  equation. 
On  the  unit  square,  we  consider 

—A u  —  k2u  =  0  on  D  =  (0, 1)  x  (0, 1)  .  . 

dnu  +  iku  =  g  on  dQ,  '  ' 

where  g  is  chosen  such  that  the  exact  solution  is  a  plane  wave  of  the  form 

7T 

u(x,  y )  —  exp{zfc(x  cos  6  +  y  sin  0)},  6  —  — . 

16 

The  following  two  types  of  local  approximation  spaces  were  analyzed  in  [11].  The 
first  type  are  “generalized  harmonic  polynomials”  as  alluded  to  in  the  introduction. 
Written  in  polar  coordinates,  they  take  the  form 

Vv(p)  =  span {e±in8Jn(kr)  \  n  =  0, . . .  ,p},  (21) 
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where  the  functions  Jn  are  the  Bessel  functions  of  the  first  kind  (see,  e.g.,  [6]).  The 
second  type  are  systems  of  plane  waves  given  by 

W{jp)  =  span{exp(iA(scos^-  +  ysmOj)) \9j  =  yj,j  =  0,...  ,p-  1}. 

(22) 

Remark  7:  The  spaces  Vv (p),  the  spaces  of  “generalized  harmonic  polynomials”, 
share  the  optimality  properties  of  the  harmonic  polynomials  for  the  approximation 
of  harmonic  functions  on  discs  (see  section  4.2);  the  spaces  Vv{p )  are  optimal  in 
the  sense  of  n- width  for  large  classes  of  rotationally  invariant  spaces  of  solutions  of 
Helmholtz’s  equation  on  discs. 

Remark  8:  The  numerical  examples  below  are  based  on  the  spaces  W{p).  In  all 
computations  we  chose  p  to  be  of  the  form  2  + 4m,  m  €  N0,  so  that  the  exact  solution 
of  problem  (20)  is  not  an  element  of  the  PUFEM  space. 

The  approximation  properties  of  these  two  types  of  spaces  are  very  similar  to  the 
usual  harmonic  polynomials.  In  fact,  we  have 

Theorem  4  Let  Cl  C  R2  be  a  simply  connected,  bounded  Lipschitz  domain.  Let 
Cl  DD  Cl  and  assume  that  u  £  L2(Cl)  solves  the  homogeneous  Helmholtz  equation  on 
Cl.  Then 

^ 11“  -  “»IU‘(0)  <  Ce-” 

where  C ,  C ,  7,  and  7  depend  only  on  Q}  fi. 

Remark  9:  For  the  solution  of  the  model  problem  (20),  theorem  4  can  be  strength¬ 
ened: 


,  llu  -  Up||jP(n)  <  C{i,Cl)e  7P, 

Vv(p ) 

.  Ilu  “  ™rll*l(n)  <  C(rf,  fi)e_7P 

wp  ew(p) 


holds  for  any  fixed  7  >  0. 


Theorem  5  Let  Cl  be  a  simply  connected  bounded  Lipschitz  domain,  star-shaped  with 
respect  to  a  ball.  Let  the  exterior  angle  of  Cl  be  bounded  from  below  by  Air,  0  <  A  <  2. 
Assume  that  u  £  Hk(Cl),  k  >  1,  satisfies  the  homogeneous  Helmholtz  equation.  Then 


j  =  0,...  ,[fc], 
j=0,...,[k}. 
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le  1:  DOF  necessary  i 

;o  obtain  accuracy  e  in  L 2  norm;  k  = 

e 

best 

approximant 

QSFEM 

GLSFEM 

FEM 

30% 

2.045D+3 

3.969D+3 

2.016D+4 

7.784D+4 

10% 

5.041D+3 

1.000D+4 

6.150D+4 

2.352D+5 

5% 

8.464D+3 

1.960D+4 

1.274D+5 

4.692D+5 

Table  2:  DOF  necessary  to  achieve  various  accuracies  in  L2  for  PUFEM  with  n  =  4 
and  various  other  methods;  A:  =  100 


V 

L 2  error 

PUFEM 

best  approx. 

QSFEM 

FEM 

i 

6.50D+2 

3.80D+3 

7.95D+3 

2.08D+5 

m 

7.50D+2 

5.89D+4 

1.23D+5 

3.23D+6 

34 

0.11% 

8.50D+2 

3.45D+5 

7.23D+5 

1.90D+7 

The  PUFEM  can  be  based  on  either  approximation  space.  In  the  numerical  results 
below,  we  concentrate  on  the  PUFEM  based  on  the  spaces  W(p)  of  plane  waves  (for 
a  comparison  with  the  “generalized  harmonic  polynomials”  Vv(p),  see  [11]).  The 
domain  Q  is  covered  by  square  patches  and  the  partition  of  unity  consists  again  of 
piecewise  bilinear  functions  as  described  in  section  5.  The  local  approximation  spaces 
Vj  are  taken  as  the  spaces  W(p). 

Remark  10:  The  theorems  above  merely  address  the  issue  of  approximability; 
we  do  not  deal  with  the  delicate  question  of  stability  of  the  finite  element  methods 
based  on  these  spaces.  Suffice  it  to  say  that  the  spaces  created  by  the  PUFEM  are 
stable  under  the  assumption  that  the  mesh  size  h  is  sufficiently  small  with  respect  to 
the  wave  number  k  (see  [11]).  However,  as  can  be  seen  in  the  numerical  results,  the 
PUFEM  performs  very  well  as  a  p  version  for  very  coarse  meshes. 

In  tables  1-6  the  PUFEM  is  compared  with  the  usual  Galerkin  finite  element  method 
(FEM),  the  generalized  least  squares  finite  element  method  (GLSFEM)  of  [22],  and 
the  quasi-stabilized  finite  element  method  (QSFEM)  of  [19].  Since  all  three  methods 
are  based  on  piecewise  linear  functions  on  uniform  grids,  tables  1  and  2  include  the 
piecewise  linear  best  approximant  for  reference.  The  FEM,  GLSFEM,  and  QSFEM 
differ  in  their  choice  of  the  bilinear  form.  In  particular,  the  bilinear  form  of  the 
QSFEM  is  constructed  such  that  “pollution”  (see  [19])  is  minimized,  and  thus  the 

Table  3:  number  of  operations  using  band  elimination  -  the  p  version  of  the  PUFEM; 
n  =  4;  k  =  100;  error  in  L2 


V 

L 2  error 

PUFEM 

QSFEM 

FEM 

26 

10.8% 

1.76D+7 

6.3D+7 

4.3D+11 

30 

0.69% 

2.71D+7 

1.5D+10 

1.01D+13 

34 

0.11% 

3.94D+7 

5.2D+11 

3.6D+14 
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Table  4:  number  of  operations  for  h/p  version  of  PUFEM;  k  =  100;  L2  error 


V 

n 

L2  error 

NOP  PUFEM 

26 

4 

10.8% 

1.76D+7 

18 

8 

10.6% 

5.23D+7 

14 

16 

9.5% 

2.75D+8 

Table  5:  operation  count  for  solving  linear  system;  error  in  H1  norm;/s  =  32 


Galerkin 

QSFEM 

VDOF 

H 1  error 

No.  iter 

NOP 

H 1  error 

No.  iter 

NOP 

hline  32 

65% 

232 

4.51D+6 

30.5% 

272 

5.29D+6 

64 

21.7% 

434 

3.37D+7 

14.3% 

492 

3.82D+7 

128 

8.16% 

831 

2.68D+8 

7.02% 

953 

2.96D+8 

256 

3.64% 

1665 

3.48% 

1863 

2.31D+9 

512 

1.72% 

3263 

1.69% 

3752 

1.86D+10 

Table  6:  operation  count  for  band  elimination  for  PUFEM;  k  =  32,  error  in  if1;  n  =  1 


V 

H1  error 

NOP  PUFEM 

■El 

46% 

1.3D+5 

2.3D+5 

26 

0.38% 

3.8D+5 

30 

0.00025% 

5.9D+5 
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QSFEM  is  virtually  the  best  method  available  which  is  based  on  piecewise  linear 
functions.  We  will  see  that  the  PUFEM  compares  very  favorably  with  the  QSFEM. 
We  discuss  the  cases  k  =  100  with  the  L2  norm  as  the  error  measure  and  k  =  32  with 
the  H 1  as  the  error  measure.  Tables  1-4  show  the  performance  of  the  PUFEM  in 
comparison  with  the  other  methods  for  k  =  100  and  the  L2  norm  as  error  measure. 
Table  1  shows  the  number  of  DOF  needed  to  achieve  a  certain  L 2  accuracy  for  the 
various  piecewise  linear  methods.  We  see  that  the  QSFEM  needs  2  times  as  many 
DOF  as  the  best  approximant,  while  the  GLSFEM  needs  10-15  and  the  FEM  40-50 
as  many.  Table  2  shows  that  the  p  version  of  the  PUFEM  can  achieve  the  same  accu¬ 
racy  as  the  other  methods  with  markedly  fewer  DOF.  This  can  be  attributed  to  the 
exponential  approximability  of  the  PUFEM:  According  to  remark  6.2  the  approxima¬ 
tion  properties  of  the  PUFEM  space  based  on  plane  waves  W(p)  are  exponential  in  p 
whereas  the  h  versions  of  the  piecewise  linear  methods  can  only  have  algebraic  rates 
of  convergence.  This  explains  why  the  discrepancy  between  the  PUFEM  and  the 
other  methods  becomes  more  pronounced  for  high  accuracy:  in  order  to  achieve  10% 
accuracy  in  L2,  the  best  approximant  needs  6  times  as  many  DOF  as  the  PUFEM, 
whereas  it  needs  400  times  as  many  as  the  PUFEM  to  achieve  0.11%  accuracy.  Ta¬ 
ble  3  shows  how  this  reduction  of  DOF  translates  into  a  reduction  in  the  operation 
count  if  a  direct  solver  (band  elimination)  is  used.  Again,  the  PUFEM  outperforms 
the  QSFEM  and  the  FEM  for  the  case  of  10.8%  accuracy  and  is  greatly  superior  for 
high  accuracy. 

In  tables  1-3  we  saw  the  performance  of  the  PUFEM  as  a  p  version.  Table  4  shows 
the  performance  of  the  PUFEM  as  an  hp  version  by  listing  the  number  of  operations 
for  the  band  elimination  for  various  combinations  of  p  and  h  =  -  which  result  in  an 
accuracy  of  ca.  10%  in  L2.  We  see  that  the  operation  count  increases  with  n  (and  thus 
with  decreasing  p).  This  can  again  be  explained  by  the  fact  that  the  PUFEM  spaces 
feature  exponential  approximability  as  p  versions  but  only  algebraic  approximability 
as  h  versions. 

Tables  5-6  illustrate  the  case  k  =  32  with  the  H1  semi  norm  as  the  error  measure. 
The  linear  system  of  the  usual  FEM  and  the  QSFEM  is  solved  using  the  iterative 
scheme  proposed  in  [4].  We  compare  the  cost  of  these  iterative  schemes  (table  5)  with 
the  cost  of  the  band  elimination  for  the  PUFEM  (table  6)  as  a  p  version  (n  =  1). 
We  see  that  the  PUFEM  is  cheaper  than  the  QSFEM,  which  is  virtually  the  optimal 
method  for  piecewise  linear  ansatz  functions.  The  PUFEM  is  cheaper  in  the  whole 
range  of  accuracies  (50%-0%).  As  in  the  case  of  DOF  versus  L2  accuracy  above, 
the  disparity  between  the  PUFEM  and  the  other  methods  becomes  bigger  for  high 
accuracy:  for  50%  error,  the  PUFEM  is  30  times  cheaper  than  the  FEM,  and  for  1% 
the  PUFEM  is  105  times  cheaper! 

Remark  11:  In  the  operation  count  for  the  PUFEM  (tables  4,  6)  only  the  con¬ 
tributions  of  the  band  elimination  are  reported.  This  is  justified  by  the  particular 
structure  of  the  problem  under  consideration.  The  mesh  is  uniform,  the  partition  of 
unity  consists  of  piecewise  bilinear  functions  and  the  local  approximation  spaces  are 
spaces  of  plane  waves.  All  of  this  can  be  exploited  in  the  construction  of  the  stiffness 
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matrix,  and  the  resulting  cost  of  the  generation  the  stiffness  matrix  is  of  lower  order 
compared  with  the  cost  of  the  linear  solver. 

The  numerical  examples  show  that  the  PUFEM  performs  much  better  than  the  usual 
h  versions  both  in  terms  of  error  versus  DOF  and  error  versus  operation  count.  This 
is  due  to  the  fact  that  the  PUFEM  allows  us  to  use  local  approximation  spaces  that 
capture  the  local  behavior  of  the  solution  very  well,  even  if  the  solution  is  rough.  In 
this  example,  the  approximation  with  plane  waves  is  very  efficient  although  the  wave 
number  k  is  large  ( k  =  32,  k  =  100).  We  saw  that  the  PUFEM  outperforms  the  h 
version  for  accuracies  of  practical  interest  (50%-l%  in  H1,  say)  and  that  the  PUFEM 
is  immensely  superior  for  high  accuracy. 

7  A-posteriori  Error  Estimation 

A-posteriori  error  estimation  for  finite  element  solutions  obtained  by  the  PUFEM  is 
possible  if  local  problems  on  the  patches  ft*  H  ft  can  be  solved  (or  suitably  approxi¬ 
mated).  In  order  to  demonstrate  this,  let  us  consider  the  model  problem 

Lu  =  —  diva(x)gradu  +  c(x)u  =  /  £  L2(Q)  on  ft 

u  =  0  on  Tn  C  ^  0 

anu  =  a(x)dnu  =  g  £  on  VN  =  dtt\  TD  (23) 

where  a,  c  are  bounded  functions  and  satisfy  the  inequality 

0  <  a  <  min(o(s),  c(z))  <  max(a(x),  c(x))  <  /3  <  oo. 

The  weak  form  of  this  problem  is  to  find  u  £  such  that 

B(u,  v )  =  F(v)  Vu  £  J2£(fl)  =  {v  £  H\n)  |v  =  0on  TD}  (24) 

where  the  bilinear  form  B  and  the  linear  functional  F  are  defined  in  the  standard 
way.  The  conditions  on  the  coefficients  a,  c  imply  that 

all‘ullif1(n)  —  B{uiu)> 

|S(u,u)|  <  /3|M|tfi(n)IM|tfi(f>)- 

Let  VFe  be  a  conforming  PUFEM  space,  i.e.,  VFE  C  Then,  the  finite  element 

solution  ufe  €  VFe  is  defined  by 

B(ufe,v)  =  F(v)  Vu  £  VFE  C  H^Q).  (25) 

On  each  patch  fi,  fl  fl,  we  introduce  the  local  problem 

find  rji  £  Wi  B(rji,v )  =  B(u  -  uFE,  v )  Vu  £  W{  (26) 
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where 


Wi  =  {v  6  H\cit  n  fl)  |  u  =  o  on  <9(fl,  n  fl)  \  vN}.  (27) 

Remark  12:  If  we  require  the  PUFEM  space  to  be  of  degree  2,  i.e.,  if  VFE  C 
He(CI)  n  (72(fl),  integration  by  parts  allows  us  to  express  the  right  hand  side  of  (26) 
explicitly  in  terms  of  the  given  data  L,  f,  and  g: 

B(r)i,  v )  =  B(u  -  uFE,  v)=  (/  -  LuFE)vdx  + 

./n,nn 

In  this  last  integration  by  parts  argument,  we  made  use  of  the  assumption  VFF  C 
C2(Cl).  This  is  an  important  simplification  in  practice  because  in  that  way,  the 
evaluation  of  the  right  hand  side  of  (26)  requires  only  knowledge  about  uFF  and  its 
derivatives  on  the  patch  fl,Tlfl.  If  the  space  VFE  is  less  regular  (e.g.,  VFE  C  C(Cl)  and 
piecewise  C2)  the  integration  by  parts  argument  introduces  additional  terms  related 
to  the  jumps  of  derivatives;  restricting  ourselves  to  the  case  VFE  C  C2(Cl)  removes 
the  necessity  to  determine  the  points  where  these  jumps  may  occur. 


/  (g  ~  crnuFE)vds. 

Jtn  (28) 


Before  we  proceed  to  prove  theorem  6,  which  relates  the  error  of  the  finite  element 
solution  to  the  local  functions  rji,  we  need  to  impose  some  approximation  properties 
on  the  local  approximation  spaces  Vi. 

Definition  3  A  collection  Vi  of  local  approximation  spaces  has  the  uniform  Poincare 
property  if  there  is  Cp  >  0  independent  of  i  such  that 

1.  for  i  such  that  Cli  fl  TE  =  0,  Vi  contains  the  constant  functions  and 
infAeR  ||v  —  ^||z.2(n,nn)  ^  Cp  diam(f2i)||u||ffi(Qt.nn)  Vu  £  H 1(fii  fl  fi) 

2.  for  i  such  that  Qi  D  Pd  ^  0 

IMU^n.-nn)  ^  CP  diam(Oi)||u||Hi(n.nn)  Vu  £  {v  £  \  v  =  0  on  I'd} 

Theorem  6  Let  {fij}  be  a  cover  of  Cl  and  a  (M,  Coo,  Cg)  partition  of  unity  sub¬ 
ordinate  to  the  cover  {f2i}.  Let  the  local  approximation  spaces  {V^}  have  the  uniform 
Poincare  property  and  assume  that  Vi  —  0  on  TD  for  Vi  £  Vi  with  Then 

there  is  C  =  C (a,  /3,  M,  C^,  Cg,  Cp)  >  0  ( which  is  explicitly  available  from  the  proof 
below)  such  that 


<  ||tt  —  UFjsIlif^n)  <  C 


(29) 


where  uFE  and  rj,  are  defined  in  equations  (25),  (26). 


Proof:  The  proof  follows  very  closely  [1].  First  we  observe  that  the  finite  element 
space  VFE  constructed  by  the  PUFEM  is  conforming,  i.e.,  VFE  C  Further¬ 
more,  we  have  Wi  C  by  continuing  the  elements  of  Wi  by  zero  on  Cl  \  (Cl,  D 12). 
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By  the  coercivity  of  B,  the  orthogonality  relation  satisfied  by  uFE,  and  the  fact  that 
ipi  =  1  on  0,  we  have 

a||n  -  ttMs||!ri(n)  <  B(u  -  uFE,  u  -  uFE ) 

=  B(u  —  Ufe,  U  —  Ufe  ~~  vFe)  \fvpE  £  VfE 

=  B(u  —  uFE,  E  <pi(u  —  uFE  —  Vi))  vFE  =  '^F(piVi,  ViGVi 
i  i 

=  ^2B(rn,(pi(u-uFE  -  v^) 


1/2 


1/2 


—  P  I  E  H^»llH1(n  t°n)  )  E  iiw(«  -  uFE  -  v<)llfrl(ninn) 


where  we  made  use  of  the  fact  that  c fi{u  —  uFE  —  vt)  6  Wi  C  Hp(Q.).  The  uniform 
Poincare  property  gives  the  existence  of  u;  £  R  such  that 

||u  -  fi|jx,2(ninn)  <  min(l,Cpdiam(fit))||u  -  UFs|Ui(n,nn)> 

and  thus  we  can  estimate 

E  ~  -  u*)llHl(ninn)  ^  E  —  Ufe  ~  u*lli2(ninn)  + 

*  i 

+  2C£,||V(u  -  UFE  -  Vt)llla(ninn)  + 

c2 

+  2diam(ni)3ll“  "  UFE  ~ 

—  /  ! 1  +  2 CqC],)  ||u  —  ^feII jf'fn.nn) 

t 

<  M  (3 Cl  +  2CqCf)  Hti  -  uF£|&1(n) 

where  we  used  lemma  2  below.  This  gives  the  upper  estimate  of  (29).  For  the  lower 
estimate,  we  use  the  fact  that  each  rji  £  Wi  C  and  thus 

E  ll77*lltfl(ninn)  ^  a  1  E  -®(77‘>77*) 

t  i 

<  or1  B(u  —  uFE,T]i)  =  B(v.  —  ufe.  ^  "  Vi) 

E* 


< /3a  |ju  -  uFF||Hl(n) 


<  /3a  1||t*  -  WFs||fri(n)Vjif  (  E  ll^llirl(ninn) 


1/2 


where  we  made  again  use  of  lemma  2  below.  This  concludes  the  proof  of  theorem  6. 

□ 
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Lemma  2  Let  LI  be  an  open  set,  {fl;}  be  an  open  cover  of  LI  satisfying  the  pointwise 
overlap  condition 

card{i  |  x  €  Lli}  <  M  Vx  £  LI. 

Let  u,  Ui  6  if1  (ft)  be  such  that 

suppit;  C  closure(ft;  D  ft). 


y~!  IMIff*(ft»)  —  -^llull A:  —  0, 1, 

i 

II  u*llH*(ft)  —  llu*llH*i(nnni)  A:  =  0, 1. 


Proof:  We  will  show  the  case  A:  =  0,  the  case  k  =  1  being  handled  similarly.  Let  Xi 
be  the  characteristic  function  of  the  domain  Lli  fl  LI.  Then 

/  M2  =  ]C  / xM2  =  [  ^2xM2<m  f  \u\2. 

i  Jcii  nn  ;  In  Jn  t  J n 

This  proves  the  first  estimate.  For  the  second  estimate,  we  use  the  overlap  condition 
and  the  condition  on  the  supports  of  the  functions  Ui  to  see  that  for  each  x  £  LI,  the 
sum  on  the  left  hand  side  extends  over  not  more  than  M  terms.  Therefore, 


/«(?“•)  w 


L2(tynty' 


Remark  13:  The  proof  of  theorem  6  shows  that  the  uniform  Poincare  property 
could  be  weakened.  It  is  enough  that  the  L2  projections  II;  :  if1  (ft;  fl  LI)  —*■  V;  satisfy 

||«  -  II;u||L2(n.nn)  <  C'pdiam(0;)||tt||Hi(n.nn). 


Remark  14:  The  existence  of  the  uniform  Poincare  constant  is  related  to  a  certain 
uniformity  of  shapes  of  the  patches.  More  precisely,  for  any  bounded  domain  D ,  the 
constant  A,  defined  by 

r>'2=  sup  inf  f t,(P> 

ueifi(£>)^R  II  Vu||x,2(25) 

is  the  second  (i.e.,  the  first  non-zero)  eigenvalue  of  the  Neumann  problem 

on  D, 
on  dD. 


—Ait 

dnu 


Ait 

0 
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Remark  15:  Let  us  note  that  a  simple  scaling  argument  shows  that  the  uniform 
Poincare  constant  of  definition  3  depends  only  on  the  shape  of  the  patches  and  not 
the  diameters.  Thus,  one  simple  way  to  enforce  a  uniform  Poincare  property  is  to 
restrict  the  number  of  possible  shapes  of  the  patches  fl;  f~|  f 1. 

Let  us  outline  sufficient  conditions  on  the  patches  fl;  fl  fl  that  guarantee  the  uniform 
Poincare  property  of  the  local  approximation  spaces  Vt  based  on  the  following  lemma. 

Lemma  3  Let  fl  C  Mn  be  a  convex  domain,  u  6  i?1( fl).  Then 

llu  -  m||l2(D)  <  (diam(fl))n  ||Vu||i2(n)  (30) 

where  u  is  the  average  of  u  over  fl 

-=w\Lu 

|fl|  stands  for  the  volume  of  fl,  and  uin  is  the  surface  of  the  unit  sphere  in  R". 

Proof:  Section  7.8  of  [5].  □ 

For  patches  fl;  D  fl  such  that  fl;  fl  Td  =  0,  lemma  3  gives  the  uniform  Poincare 
property  if  fl;  fl  fl  is  convex  and  if  there  is  p  >  0  such  that  each  patch  contains  a  ball 
of  radius  pdiam(fl;)  (and  is  trivially  contained  in  a  ball  of  radius  diam(fl;)).  Note 
that  this  is  reasonable  restriction  on  the  patches  in  view  of  condition  (4). 

Let  us  now  turn  to  the  patches  close  to  the  boundary  where  the  Dirichlet  conditions 
are  prescribed.  For  simplicity,  consider  a  two  dimensional  setting,  assume  that  the 
patches  fl;  are  discs,  and  that  fi;  fl  rD  ^  0,  fl;  fl  =  0.  Moreover,  let  fl;  fl  TD  be  a 
straight  line  segment.  If  fl;  fl  fl  is  less  than  a  half-disc  (but  fl,  n  fl  still  contains  a  ball 
with  diameter  /5diam(fl,)),  the  reflection  across  T#  yields  a  convex  domain  fl  C  fl*. 
For  u  fl  fl)  such  that  u|rD  =  0,  the  antisymmetric  extension  across  To  gives 

an  H1( fl)  function  with  zero  average,  and  thus  lemma  3  gives 

\/2 

2|M|z,2(n<nn)  =  ||w||x,2(n)  <  —  diam(fli)||Vu||z,2(ninn)- 

The  case  that  fl*  fl  fl  is  bigger  than  a  half-disc  can  be  reduced  to  the  above  one  by 
an  appropriate  mapping.  The  necessary  condition  for  the  Poincare  constant  not  to 
degenerate  is  that  the  length  of  the  line  segment  fl;  fl  Tp  >  pdiam(fl;). 

The  case  fli  Pi  Tx?  ^  0,  fl;  fl  Tjv  yf  0  can  be  dealt  with  using  similar  ideas.  Again,  the 
necessary  condition  is  that  the  length  of  the  line  segment  fl;  fl  >  pdiam(f2;). 
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